%%%
%% Sorry, TEX, I'll be back! I promise to rewrite the thing after finall 
%% approvement of the texts, temporarily written using MSWord. February 5th, 2011
%%
%% As I've promised ^^. February 24th, 2011
%%%
\documentclass[a4paper,12pt]{article}
\usepackage{mathtext}%Кирилица в математическом режиме
\usepackage[russian]{babel}
\usepackage[utf8]{inputenc}
\usepackage{amssymb}
\usepackage{amsmath}%including cases
\usepackage{graphicx}
\renewcommand{\baselinestretch}{1.2}
\usepackage[top=2cm,text={160mm,255mm},centering]{geometry}


%\textwidth=160mm \textheight=255mm 
%\oddsidemargin=5mm
%\topmargin=-10mm
\parindent=1cm
%\parindent=0em
%\mathsurround=2pt \pagestyle{empty} \sloppy

%Макросы для удобства
%\newcommand{\rot}{\mathbf{rot}}
%\newcommand{\diverg}{\mathbf{div}}
%\newcommand{\ksi}{\xi}
%\newcommand{\citepage}[2]{\cite[стр. #1]{#2}}
%\newcommand{\eps}{\varepsilon}

%\newcommand{\inttraj}{траект.}%{\sim}
%\newcommand{\intt}{\int\limits_{\inttraj}dl ~~}
%\newcommand{\refp}[1]{(\ref{#1}) на стр. \pageref{#1}} % Ссылка далеко
%\newcommand{\refl}[1]{(\ref{#1})} % Ссылка близко
%\newcommand{\leqR}{\leqslant}%Russian \leq
%\newcommand{\geqR}{\geqslant}%Russian \geq
%\newcommand{\atval}[2]{\bigg|_{#1}^{#2} }

%\newcommand{\pdi}{\partial_i}
%\newcommand{\pdj}{\partial_j}
%\newcommand{\pdk}{\partial_k}
%\newcommand{\pdeta}{\partial_\eta}





\newcommand{\parsec}{\hspace{1cm}}


\newcommand{\ilus}[3]{
	\begin{figure}[h!]
		\centering
		\includegraphics[width=15cm]{#1}
		\caption{\footnotesize\bfseries #2}
		\label{#3}
	\end{figure}
	{\itshape\bfseries См. график под номером (\ref{#3}) на странице (\pageref{#3}).}
}
%\newcommand{\ilus}[2]{
	%\begin{center}
		%\includegraphics[width=13cm]{#1}\\
		%{\footnotesize\bfseries #2}
	%\end{center}
%}

\begin{document}
{
\pagestyle{empty} % отключает нумерацию
\large
	\begin{center}
	\vspace*{10mm}
	Санкт-Петербургский Государственный Университет\\
	Физический факультет\\
	Кафедра радиофизики\\
	\vspace{49mm}
	Магистерская диссертация\\
	Студента 7 курса\\
	Зелинского Ивана Сергеевича\\
	«Электромагнитная модель молниевого разряда.»\\
	\vspace{45mm}
	Научный руководитель: доцент, кандидат физмат наук,\\
	Кононов Игорь Иванович.\\
	\vspace{64mm}
	{
		\small
		Санкт-Петербург\\
		2011 г.\\
	}
	\end{center}
	\newpage
}
\tableofcontents
\newpage

{
	\centering
	\large
	\bfseries
	Электромагнитная модель моль молниевого разряда\\
}

\section{Введение}
\parsec Описание и моделирование разрядных процессов, возникающих при возникновении  молниевой вспышки, а также сопровождающего ее развитие электромагнитного излучения, является объектом исследований, проводимых с момента открытия молнии как электрического процесса. В свою очередь, из очень сложного по своей пространственно-временной структуре импульсного излучения молниевых вспышек, охватывающего широкий частотный спектр от долей Гц до сотен МГц, особый интерес вызывают их сильноточные импульсные компоненты. К ним можно отнести толчкообразно (ступенчато) развивающиеся лидерные процессы, характеризуемые пространственными размерами от десятков до сотен метров и  токами, достигающими единиц кА, первые и повторные обратные удары (иногда называемые главными разрядами). Каждый такой разряд начинается с развития направленного вверх относительно короткого (несколько десятков метров) стримера при приближении стримерной зоны  головки спускаемого лидера к поверхности земли. После их смыкания начинается фаза мощного импульсного пробоя в виде ярко светящегося разряда, направленного к грозовому облаку. Токи наиболее сильноточных обратных ударов, протекающие по плазменным каналам, образованным предшествующим лидерным процессом (с диаметром до нескольких метров),  достигают десятков килоампер. Стягиваясь под действием собственного магнитного поля в очень узкий плазменный шнур  с диаметром, составляющим  несколько сантиметров. В начале своего развития обратные удары имеют преимущественно вертикальную  ориентацию, переходящую по прошествии 30-40 микросекунд в протяженную внутриоблачную часть преимущественно горизонтальной ориентации. У некоторых специфических разрядов (называемых «положительными» по знаку заряда, переносимого к земле) пробой начинается с положительно заряженной части облака, располагающейся на несколько км выше области отрицательного заряда, из которой обычно стартует ступенчатый лидерный процесс. Более значительный разрядный промежуток между облаком и землей у «положительных разрядов» приводит к более значительным разрядным токам. К тому же для них характерна более протяженная горизонтальная составляющая. В силу этих причин такие разряды считают (пока гипотетически) стимуляторами возникновения джетов и спрайтов, развивающихся в стратосфере и нижней части ионосферы. 

К классу сильноточных разрядов можно отнести  также внутри- и межоблачные стримеры, развивающиеся между пространственно разделенными положительно и отрицательно заряженными областями, суммарный заряд которых достигает нескольких единиц и даже десятков кулон. Их развитие сопровождается появлением последовательности (цуга)  так называемых К-импульсов, амплитудные распределения которых перекрываются с нижней частью распределений амплитуд обратных ударов между облаком и землей.

Все перечисленные выше типы разрядов имеют спектральную плотность радиочастотной составляющей, сосредоточенную в ИНЧ-ОНЧ диапазонах с максимумами, лежащими в интервале 3-5 кГц для обратных ударов, около 20 кГц для К-импульсов и 30-50 кГц для лидеров. Амплитудные значения импульсов напряженности, например, электрического поля достигают на удалении 100 км соответственно 5-10В/м у обратных ударов и 1-2В/м у К-импульсов и лидеров. Эти особенности обеспечивают распространение этих импульсов на очень большие расстояния с относительно малым затуханием, обеспечивающим  их прием с высоким отношением сигнал-шум вплоть до расстояний в несколько тысяч километров. Все это обуславливает большую практическую значимость этого типа источников импульсного электромагнитного излучения для решения многих прикладных задач и, как следствие, непрекращающийся интерес к разработке адекватной модели сильноточного молниевого разряда как источника этого излучения. Разработка такого рода модели, по возможности наиболее полно учитывающей пространственно-временные особенности реальных сильноточных молниевых разрядов составляет цель и основное содержание настоящей работы.


\section{Обзор существующих моделей.}
\parsec В настоящее время существует большое число публикаций (превышающее сотню наименований), посвященных разработке моделей молниевых разрядов. Приведем краткий обзор существующих моделей, делая акцент на тех из них, которые наилучшим образом характеризуют их в качестве источников ЭМИ, позволяющих с наибольшей простотой и в то же время с максимально достижимой точностью использовать при проведении расчетов полей на различных удалениях от излучателей в меняющихся условиях распространения.

Следуя определениям работы [1], можно выделить несколько основных классов (типов) моделей сильноточных молниевых разрядов, рассматриваемых и обсуждаемых в научной литературе.

Первый класс так называемых газодинамических или «физических» моделей связывают с  радиальной эволюцией короткого сегмента молниевого канала и сопровождающей динамику его развития ударной волной, описываемых системой газодинамических уравнений, отражающих основные законы сохранения (массы, момента количества движения, энергии). Их решение позволяет оценить температуру, давление и массовую плотность как функцию времени и радиальной координаты. Некоторые модели этого типа [3-5] были развиты в приложении к лабораторным искровым разрядам, но без достаточных на то оснований применяются для обратных ударов молниевого происхождения.

Второй класс моделей, называемых электромагнитными основываются на представлении разряда в виде тонкой проволочной антенны с потерями. Возбуждаемый в антенне ток рассчитывается численным решением уравнений Максвелла. В модели [7] рассматривается прямой вертикальный канал, в то время как авторы модели [6] имеют дело с пространственно распределенным каналом произвольной формы, включающим ветвления, поражаемые объекты и учитывают нелинейные эффекты в «соединительных» процессах. Автор работы [8], используя уравнения Максвелла, описал процессы в стреловидном лидере и обратном ударе как направленные волны, распространяющиеся вдоль проводящего цилиндрического канала. 

Третий класс  моделей является, по сути, аппроксимацией модели предыдущего класса в виде вертикальной длинной линии (transmission-line model) с распределенными параметрами. Пространственно-временное распределение тока вдоль канала рассчитывается как переходной процесс, возбуждаемый  заданным током в его начале.

Параметры импульсов электромагнитных полей (величина начального пика, времена нарастания и спада, значение нулевого перехода), рассчитанные в работах [9] в рамках линейно распределенной модели и с учетом нелинейных эффектов в работе [10], сопоставлялись с аналогичными характеристиками экспериментальных форм, приведенными в [11], не дали удовлетворительных результатов.

Не обсуждая далее детали перечисленных выше первых трех классов моделей, отметим лишь, что их сравнительно каткий обзор и анализ приведены в работе [1].

Четвертый класс так называемых «инженерных» моделей, начало рассмотрения которых было положено работой [2], предполагает задание пространственно-временного распределения тока на основе наблюдаемых (измеряемых) типичных форм токов в основании молниевого канала, профиля его светимости и скоростей движения фронта. В этом классе моделей физическая сторона явлений сознательно не рассматривается, основной упор делается на подбор формы и параметров переходного тока, по- возможности, адекватно описывающих рассчитанные с его использованием электромагнитные поля на расстояниях от нескольких десятков метров до сотен километров от основания разряда. Подбор формы тока и уточнение его параметров осуществляется путем сопоставления расчетных полей с экспериментально наблюдаемыми.

Методика подтверждения применимости рассматриваемых моделей основывалась на двух несколько отличающихся подходах. В одном из них для расчетов и последующего сопоставления использовалась некоторая типовая усредненная форма тока в основании разряда с усредненными же значениями описывающих ее параметров, полученных в ряде экспериментальных регистраций при прямых ударах молний в горных обсерваториях или на высотных объектах. Расчетные формы затем сравнивались с типовыми, также усредненными формами импульсов, зарегистрированными на различных расстояниях от источников излучения. 

В другом, более точном подходе, использовалась конкретная форма тока в начале разряда, полученная от молнии, инициируемой выносимой вверх с помощью небольшой ракеты металлической проволокой. Разряд происходил по достижении потенциала приближающегося грозового облака некоторого критического значения, когда проволока раскручивалась до высоты в несколько сот метров. По полученной форме тока с использованием скорости продвижения ее фронта, оцениваемой с помощью аппаратуры оптической регистрации, рассчитывались поля на расстояниях, соответствующих позициям приемников, регистрирующих ЭМИ синхронно с развитием инициируемого разряда. Здесь удавалось подобрать параметры волны тока, распространяющейся вдоль канала, с минимальным расхождением расчетных полей с экспериментально зарегистрированными. Основным  недостатком этой методики является то обстоятельство, что формы и параметры разрядных токов и соответствующих им полей для триггерных молний существенно отличаются от молний и полей обычных равнинных гроз.

Характерной особенностью класса «инженерных» моделей является относительно небольшое число оцениваемых параметров, что определило его широкое распространение  и обсуждение в литературе. Таким образом, инженерные модели, при всей своей упрощенности, являются в то же время наиболее согласованными с измеряемыми электрическими и магнитными полями от естественных и инициируемых (триггерных) молний. 

К настоящему времени разработано большое число различных версий инженерных моделей, отличающихся зачастую малосущественными добавками, призванными уточнить отдельные достаточно мелкие детали наблюдаемых форм. Не проводя обсуждения и сравнительного анализа различных подходов при разработке этих моделей, достаточно подробный обзор которых приведен в работах [1,12], мы попытаемся разработать и программно реализовать более универсальную модель инженерного класса, отказавшись от многих ограничений, вводимых в известных нам публикациях. Пользуясь этой моделью, мы проведем оценки границ применимости более простых представлений разряда, в частности, в виде точечного диполя (вертикального для обратных ударов облако-земля и произвольно ориентированного для внутриоблачных разрядных процессов). Это позволит существенно упростить методику расчетов импульсов электромагнитного излучения (ЭМИ) сильноточных молниевых разрядов в сложных условиях распространения, используя возможности теоретического и вычислительного аппарата,  разработанного на кафедре радиофизики СПбГУ.


\section{Протяженный одномерный источник в качестве модели молниевого разряда.}
\parsec В подавляющем большинстве публикаций молниевый канал рассматривается прямолинейным и вертикально ориентированным, а точка наблюдения расположена на поверхности земли. Отказавшись от этих допущений, будем рассматривать молниевый канал в виде нитевидной структуры произвольной конфигурации и протяженности, расположенной над поверхностью земли (в обычной практике моделирования рассматривается вертикальный линейный канал), а точка наблюдения может быть расположена в любом произвольном месте над поверхностью земли (обычно поля рассчитываются на ее поверхности). Это позволит рассчитывать поля практически всех видов разрядов, происходящих как между облаком и землей (лидеров и обратных ударов), так и внутри облаков (для произвольно ориентированных стримеров, генерирующих К-импульсы).

Чтобы выделить в возможно более «чистом» виде особенности и параметры ЭМИ разряда, связанные с геометрической структурой молниевого канала и параметрами распространяющегося вдоль него тока, расчет различных компонент поля будем проводить на расстояниях до 20-30 км от основания  разряда. С учетом рассматриваемого частотного диапазона (до 100кгц) Землю можно считать плоской и бесконечно проводящей.

В рамках сделанных допущений молниевый разряд описывается следующим образом:
\begin{description}
	\item[а)] геометрическая структура канала описывается кривой:
	$\overrightarrow{{{R}_{T}}}(t):[0,1]\mapsto {{\mathbb{R}}^{3}}$ длина которой $L : [0,1] \mapsto [0,L_{max}]$
	\item[б)] ток вдоль канала $I = I(t,l), t \in [0,+\infty), l \in [0,L_{max}]$, 
	определяемый в виде импульса заданной формы с изменяющейся (уменьшающейся) амплитудой, включаемый в момент времени t=0 и  распространяющийся  вдоль него с переменной скоростью, может быть записан в виде: 
	\begin{equation}I(t,l)=\Phi \left( t-\tau (l) \right)\xi (l) \label{ref_1}\end{equation}
\end{description}
Величины, входящие в формулу (\ref{ref_1}), подчиняются следующим требованиям:
	\begin{equation}\Phi(t) = 0, ~ \forall t < 0 ~~ \Phi(+\infty) = 0\label{ref_2}\end{equation}
	\begin{equation}\tau(0) = 0~~\frac{d\tau\left(l\right)}{dl} > 0\label{ref_3}\end{equation}
	\begin{equation}\xi(0) = 1~~\xi(l) > 0~~\frac{d\xi(l)}{dl} \leqslant 0\label{ref_4}\end{equation}
В большинстве публикаций отыскание поля источника подобного излучателя производится во временной области. Однако, такие методы неприменимы к задачам распространения, теория которых основывается на спектральных представлениях. В настоящей работе мы используем спектральные методы. Это позволяет при необходимости развить теорию так, чтобы воспользоваться существующими методами решения задач излучения точечного диполя в сложной среде. \\
\\
Будем искать решение, используя преобразование Фурье в виде:
\[S_w = \frac{1}{2\pi}\int\limits_{-\infty}^{+\infty}S(t)e^{i\omega t}dt\]
Разложим ток $I(t,l)$ в Фурье-спектр по времени:
\[{{I}_{w}}(l)=\frac{1}{2\pi }\int\limits_{-\infty }^{+\infty }{\Phi \left( t-\tau (l) \right)\xi (l){{e}^{i\omega t}}dt}\]
После упрощения получается:
\begin{equation}
I_w(l) = \xi(l) e^{i\omega\tau(l)}\Phi_\omega \label{ref_5}
\end{equation}
где $\Phi_\omega$ - спектр $\Phi(\tau)$.
Для конкретных численных оценок мы воспользуемся некоторыми модельными аппроксимациями, наиболее часто встречающимися в литературе, а именно:
	\begin{equation}\Phi(t) = I_0 \left( \exp(-\alpha t) -\exp(-\beta t)\right)\theta(t)\label{ref_6}\end{equation}
или
	\begin{equation}\Phi(t) = I_0~ t^n \exp(-\alpha t)\theta(t)	\label{ref_7}\end{equation}
Зададим запаздывание следующим дифференциальным уравнением:
	\[\frac{dl}{d\tau}  =  v_0 \exp(-\gamma \tau)\]
из которого следует
	\[\tau(l) =- \frac{1}{\gamma} \ln \left(1 - l \frac{\gamma}{v_0} \right)\]
Определим затухание:
	\[\xi(l) = \exp(-\eta l)\]
Здесь $\tau$ - запаздывающее время. Параметр $\eta$ ответственен за затухание волны тока, параметр $\gamma$ - за её замедление, наблюдаемое на практике в оптических измерениях светимости канала. $\theta(t)$ - функция включения. Приведённые выражения удовлетворяют условиям (\ref{ref_2}), (\ref{ref_3}), (\ref{ref_4}).

Фурье образы приведенных аппроксимаций, используемые при дальнейших расчетах, имеют простой аналитический вид, для (\ref{ref_6}):
\begin{equation}\Phi_\omega = \frac{I_0(\beta - \alpha)}{2\pi(\beta-i\omega)(\alpha-i\omega)} \label{ref_8}\end{equation}
Для (\ref{ref_7}):
\begin{equation}\Phi_\omega = \frac{I_0 n!}{2\pi(\alpha-i\omega)^{n+1}} \label{ref_9}\end{equation}
Вместо величины $I_0$ удобно использовать заряд, перенесённый током $\Phi(t)$ за время развития процесса: 
	\[Q = \int\limits_0^{+\infty} \Phi(t) dt = \Phi_\omega(0)\]
Спектр $\Phi_w$ известен: выражения (\ref{ref_8}), (\ref{ref_9}) при $\omega = 0$ связывают $Q$ и $I_0$. В результате,  (\ref{ref_8}) переписывается как:
	\[\Phi_\omega = \frac{\alpha\beta Q}{2\pi(\beta-i\omega)(\alpha-i\omega)}\]
и (\ref{ref_9}) как:
	\[\Phi_\omega = \frac{Q\alpha^{n+1}}{2\pi(\alpha-i\omega)^{n+1}}\]
Заряд более удобен для нормировки чем ток, тем более что $I_0$ - является не более чем константой размерности тока: связь её с токовым максимумом зависит от типа источника.

\subsection{О сохранении заряда.}
\parsec Применим к плотности тока $I(t,l)$ уравнение сохранения заряда.
	\[\partial_l I(t,l) + \partial_t \rho(t,l) = 0\]
Здесь $\rho$ - линейная плотность заряда, величина размерности $\left[ Q/L \right]$.
Пусть до начала процессов плотность заряда была равна нулю всюду. При каждом отдельно взятом $l$ проинтегрируем уравнение по времени:
	\[\rho (t,l)=-\int\limits_{0}^{t}{{{\partial }_{l}}}I({t}',l)d{t}'\]
Подставим обобщенное выражение (\ref{ref_1}):
	\[\rho (t,l)=-\int\limits_{0}^{t}{\left[ {\xi }'(l)\Phi ({t}'-\tau (l))-{\Phi }'({t}'-\tau (l)){\tau }'(l)\xi (l) \right]}d{t}'\]
	\[\rho (t,l)=-{\xi }'(l)\int\limits_{0}^{t}{\Phi }({t}'-\tau (l))d{t}'+\xi (l){\tau }'(l)\int\limits_{0}^{t}{{{\Phi }'}}({t}'-\tau (l))d{t}'\]
В пределе $t \rightarrow +\infty$:
	\[\rho (+\infty ,l)=-{\xi }'(l)\int\limits_{0}^{+\infty }{\Phi }({t}'-\tau (l))d{t}'+\xi (l){\tau }'(l)\int\limits_{0}^{+\infty }{{{\Phi }'}}({t}'-\tau (l))d{t}'\]
Используем свойство $\tau(l) > 0$ для упрощения пределов интегрирования и перепишем второй интеграл, взятый от производной:
	\[\rho (+\infty ,l)=-{\xi }'(l)\left( \int\limits_{0}^{+\infty }{\Phi }({t}')d{t}' \right)+\xi (l){\tau }'(l)\left( \Phi ({t}')|_{0}^{+\infty } \right)\]
Второе слагаемое равно нулю из-за условия (\ref{ref_2}).
	\begin{equation}\rho (+\infty ,l)=-{\xi }'(l)\int\limits_{0}^{+\infty }{\Phi }({t}')d{t}'	\label{ref_10}\end{equation}
Полученное уравнение позволяет определить, оставляет ли после себя модельный импульс тока заряды на траектории. 
Если форма тока такова, что переносится отличный от нуля заряд (интеграл в формуле (\ref{ref_10}) не равен нулю), только незатухающая волна ($\xi = const$) не оставляет после себя зарядов.

Стоит заметить, что вопросы, связанные с тем, остаются заряды на траектории или нет, имеют для нас второстепенное значение. В самом деле, мы не претендуем на то, чтобы описывать процессы, происходящие в канале после разряда – для нас важен мощный импульс, излучаемый в процессе протекания волны тока.  


\section{Расчет полей в ближней зоне.}
В соответствии с принятой выше  моделью для решения нашей задачи достаточно вычислить сумму полей исследуемого источника и его отражения от земной поверхности. Отраженный источник вида (\ref{ref_1}) имеет траекторию, зеркальную исходному, и ток $\Phi(t)$ противоположного знака.
Разобьём траекторию на сегменты и сосчитаем для каждого вклад в потенциал Герца: 
	\[\delta\vec{\Pi}_{e\omega} = \frac{i}{\omega} I_\omega(l)\delta l \frac{1}{4\pi} \frac{e^{i k R}}{R}\]
	\begin{equation}\vec{\Pi}_{e\omega} = \frac{i}{4\pi \omega}\int\limits_{\sim\sim} dl I_\omega(l) \vec{\tau}(l) \frac{e^{i k R}}{R} \label{ref_11}\end{equation}
Имея потенциал Герца, можно вычислить поля $\vec{E}$ и $\vec{H}$.
	\begin{equation}\vec{E}_\omega = \frac{1}{\varepsilon}\left(k^2\vec{\Pi}_{e \omega} + \nabla \mathbf{div} \vec{\Pi}_{e \omega}\right) \label{ref_12}\end{equation}
	\begin{equation}\vec{H}_\omega = -i\omega\mathbf{rot}\vec{\Pi}_{e\omega} \label{ref_13}\end{equation}
После применения дифференциальных операторов к интегралу для потенциала Герца получается:
{\small
	\begin{equation}E_{wj}(\vec{r}) = \frac{ik}{4\pi\varepsilon\omega} \int\limits_{\sim\sim}dl'
I_\omega\left(l'\right)\left[ \frac{k e^{ikR}}{R}\left(\tau_j - n_j\left(\vec{\tau},\vec{n}\right)\right) + 
\frac{ie^{ikR}}{R^2} \left(\tau_j - 3 n_j \left(\vec{\tau},\vec{n}\right)\right)\left(1 + \frac{i}{kR}\right)
 \right] \label{ref_14}\end{equation}
	\begin{equation}H_{wj}(\vec{r}) = \frac{1}{4\pi}\int\limits_{\sim\sim} dl' I_\omega\left(l'\right)\left[\vec{n},
\vec{\tau}\right]_j\frac{e^{ikR}}{R}\left(ik-\frac{1}{R}\right) \label{ref_15}\end{equation}
}
В формулах применяются обозначения:
\begin{equation}\vec{n} = \frac{\vec{R}}{R},\hspace{1cm}R = \left|\vec{R}\right|,\hspace{1cm}\vec{R}=\vec{r}-\overrightarrow{{{R}_{\text{}}}}({l}')\label{ref_16}\end{equation}
Подробности аналитических преобразований при получении приведенных выражений вынесены в {\itshape Приложение 1}.

Рассмотрим теперь некоторые детали применяемых численных методов. Для вычисления сигналов во временной области разумно использовать развитый аппарат быстрого преобразования Фурье. Это предполагает периодическое продолжение сигнала с периодом $T$, при котором возможно появление ошибок, связанных с наложением копий сигналов друг на друга. Величина ошибки зависит от поведения сигнала при времени, стремящемся к бесконечности. Свойства затухания влияют на выбор периода повторения. В используемых алгоритмах эта длительность выбирается по тройной сумме: времени прохождения импульса вдоль канала, запаздывания от максимально удаленной точки канала до точки наблюдения и периода такого, что большая часть относительного остатка интеграла от квадрата тока менее заданного малого параметра. По аналогичной методике оценивается предельная частота быстрого преобразования Фурье $\Omega_{max}$.

 $T$ выбрано так, что «хвосты» сигнала с достаточной точностью можно считать не перекрывающимися. К периодическому сигналу применяется теория рядов Фурье.
 
Рассмотрим некоторые особенности при вычислении поля $E$, связанные с наличием полюса в спектре в точке $\omega = 0$. Поле $E$, в общем, содержит в себе статический скачек. Величина поля после окончания процесса отлична от той, что была до его начала. Этому соответствует полюс $\omega = 0$ в формуле (\ref{ref_14}). Сигнал, состоящий из компонент $E$, нельзя периодизировать, как это было описано выше, но можно размножить во времени его производную - она не содержит скачка. После интегрирования полученного сигнала, поле $E$ будет восстановлено с точностью до константы. Эта неопределенность означает, что мы не можем знать статическое поле до начала процесса: в самом деле, наша модель имеет дело с токами, но не с изначальной конфигурацией зарядов. 

Можно доказать, что если приравнять коэффициент с номером $k=0$ в ряде Фурье для компонент поля $E$ (который нельзя вычислить из-за особенности), получившийся сигнал будет отличаться от искомого поля на произвольную прямую $Ct+D$. В самом деле, будем вычислять поле$E$, используя его производную, как описано выше. $E'_w = -i \omega E_w$ по теореме о спектре производной. Для получившейся функции напишем вычисляющий ряд Фурье, его коэффициенты: $c_k = -i\omega E_w (k \Omega)$.

Проинтегрируем ряд по времени для восстановления $E$, изъяв из ряда нулевой коэффициент. Подынтегральное выражение отличается от $E'_t$ на константу $C$, значит$E_t$ будет отличаться от результата не на $D$ а на $t+D$. В правой же части, после интегрирования по времени ряда Фурье, получится новый ряд с коэффициентами $c_k = E_w(k \Omega)$ при $k$ отличном от нуля и $c_k = 0$ при $k=0$. Это равенство - то, что требуется для доказательства.

Как было выше сказано, ряды Фурье, полученные из исследуемого сигнала, считаются при помощи быстрого преобразования Фурье. Быстрое преобразование Фурье размерности $n$: ${{x}_{j}}={\mbox{fft}_{n}}{{({{z}_{k}})}_{j}}$ - определяется как сумма $x_j = \sum\limits_{k=0}^{n-1}z_k e^{-2\pi i \frac{jk}{n}}$. Приведём практически полезную формулу, в которой все преобразования спектра сигнала сведены воедино. Пусть $T$ - длительность такая, что перекрытием сигналов при периодическом продолжении можно пренебречь. Пусть ${{\Omega }_{\max }}$ - такая круговая частота, что значения спектра сигнала при $\omega >{{\Omega }_{\max }}$ пренебрежимы. Тогда значения сигнала во времени восстанавливаются при помощи следующих выражений:

\[S(t_j) = \Omega ~\mbox{fft}_n~(z_k)\]
\[t_j = \frac{T}{n}j\]
\[
{z}_{k}=\begin{cases}
& k\ge 0,~k\le \mu :S(k\Omega ) \\ 
& k\ge \mu +1,~k\le n-1:S((k-n)\Omega ) \\ 
\end{cases}
\]
\[\mu = n / 2~~(без~остатка)\]
\[\Omega = \frac{2\pi}{T}\]
\[n \geqslant \frac{2\Omega_{max}}{\Omega}\]
При вычислении коэффициентов $z_k$ следует пользоваться свойством $S(-w)={{S}^{*}}(w)$ спектра вещественного сигнала - скопировать сопряженные значения уже вычисленных гармоник. Размер быстрого Фурье-преобразования $n$ следует выбирать нечётным, в этом случае величина $\mu$ разделит спектр на две равные доли, состоящие из отрицательных и положительных частот, и одной нулевой частоты. Вообще говоря, поскольку коэффициенты $z_k$ при $k$ близких к $n/2$ соответствуют высоким частотам, значение их несущественно. Это означает, что выбор границы $\mu$, немного отличный от описанного здесь, не должен изменить получившийся сигнал существенным образом.


\section{Некоторые результаты численных расчётов.}
\parsec Приведём конкретные результаты расчётов, выполненных для некоторых простейших моделей разряда, которые имеют как самостоятельный интерес, так и могут быть использованы для сопоставления с расчётами, выполненными для той же модели другими методами.   
\subsection{Простой вертикальный канал.}
\parsec Вычислим поля источника, представляющего из себя вертикальный канал высотой 3 км.  Форма тока – две экспоненты (\ref{ref_6}). Наблюдатель находится на расстоянии 10 км от точки (0,0,0), совпадающей с началом канала и находящейся на поверхности Земли. Цилиндрическая координата $\theta = \pi/2 - \alpha$ наблюдателя меняется от нуля до $\pi/2$, при этом наблюдатель не покидает плоскости $\phi = 0$. Для полей $E$ и $H$ вычисляются величины  $-E_{\theta}$ и $H_\phi$,  минус у электрического поля введён для удобства.\\
Параметры токовой формы: $\alpha = 10^4$[Гц], $\beta = 10^5$[Гц].  Суммарный заряд $Q = 3.87$[Кл]. Скорость волны неизменна и равна 0.7 скорости света, волна не затухает.
\ilus{MultiVertical/Hlinear.eps}{Поле H для простого вертикального источника.}{ris_1}
\ilus{MultiVertical/Elinear.eps}{Поле E для простого вертикального источника.}{ris_2}
Обратим внимание на излом, расположенный после основного максимума поля H и минимума поля E. Он возникает из-за различия задержки распространения от реального и мнимого источников.

Приведём отдельно формы вкладов от источника и его отражения для угла $\alpha ={{41}^{{}^\circ }}$.  Зелёная кривая – поле реального источника, синяя – отражения, красная – суммарный сигнал. 
\ilus{updown/updownH.eps}{Вклад отражения, поле H.}{ris_3}
\ilus{updown/updownE.eps}{Вклад отражения, поле E.}{ris_4}
Приведённые графики показывают различия в полях реального источника и его отражения: сигналы отличаются в два раза по амплитуде и имеют разные времена достижения максимума, что связано с задержкой распространения.

Поля, вычисленные для вертикального канала, были сопоставлены с источником [14], при этом результаты совпали. Совпадение является подтверждением правильности цепочки вычислений от теории, до практической реализации.

Изменим источник. Добавим замедление скорости продвижения импульса.
Параметр $\gamma = 3.56~ ~10^4$[Гц], угол $\alpha ={{41}^{{}^\circ }}$.  Красная кривая – сигнал без замедления, зелёная – с замедлением.
\ilus{speedfade/alpha_41deg_H_fade.eps}{Влияние замедления скорости импульса на форму сигнала, поле H.}{ris_5}
\ilus{speedfade/alpha_41deg_E_fade.eps}{Влияние замедления скорости импульса на форму сигнала, поле E.}{ris_6}
Как видно, небольшие изменения затронули только верхушку импульса: особенности проявляются при протекании быстро меняющегося фронта токовой волны в начале процесса. При этом, выбранное значение параметра $\gamma$ соответствует замедлению скорости до половины от изначального значения.

Сравним импульс, наблюдаемый на достаточном удалении, с импульсом точечного диполя. Для этого, построим полярную диаграмму величины $\int\limits_0^{+\infty} |H_\phi|^2 dt$.  Эта величина пропорциональна энергии в том случае, когда распространение носит волновой характер. Для сравнения, нормируем величину на значение её при $\theta = \pi / 2$. Полярная кривая для точечного диполя известна: $r = (sin\theta)^2$.
Расчёты приведены для вертикального канала длиной 4 км. Начальная скорость $v_0 = 0.7c$, $\gamma = 2.67~~10^4$ Остальные параметры не изменены. 
\ilus{compdipole/30km_v0.7_gamma_2.67e4.eps}{Диаграмма распределения энергии сложного источника в сравнении с точечным диполем.}{ris_7}
Приведённый график не позволяет отличить теоретическую дипольную кривую (синяя) от расчётной (красная). Ситуация не меняется при увеличении начальной скорости $v_0$ до 0.95c (приводить два одинаковых графика не имеет смысла). Полученный результат означает, что отличия в распределении энергии из-за близости скорости распространения импульса тока к скорости света не различимы при выбранной длительности токового импульса (параметры $\alpha$ и $\beta$ равные$10^4$ и $10^5$ соответственно).
\subsection{Источник сложной формы}
\parsec Получим поля для источника более сложной геометрической формы. Пусть канал состоит из трёх сегментов, задаваемых координатами точек соединения: (0,0,0), (0,0,1), (0,0.7,1.7), (0,0.7,2.7). Одна единица, используемая в этом обозначения, равняется двум километрам: для получения истинных значений следует умножить каждый из четырёх векторов на 2 км. Безразмерные единицы удобно использовать для краткости записи координат и возможности масштабирования геометрической структуры. Наблюдатель удален он начала координат на 30 км. Параметр замедления $\gamma = 1.79~~10^4$.
\ilus{zagagulja/Ezagagulja.eps}{Поле E источника сложной формы.}{ris_8}
\ilus{zagagulja/Hzagagulja.eps}{Поле H источника сложной формы.}{ris_9}
Как и следовало ожидать, получился сигнал сложной формы.

Траектория, состоящая из трех прямолинейных сегментов, является модельным приближением. Заменим кусочно-линейую траекторию гладкой кривой Безье третьего порядка. 
\ilus{zagagulja/trajectories.eps}{Кусочно-линейная и гладкая формы канала.}{ris_10}
Кривая Безье порядка $n$ задаётся выражением
\[\vec{R}(t)=\sum\limits_{k=0}^{n}{C_{k}^{n}{{{\vec{P}}}_{k}}{{t}^{k}}{{(1-t)}^{n-k}}}\]
Параметр $t$ изменяется в пределах от нуля до единицы. Опорные векторы ${{\vec{P}}_{k}}$ определяют форму кривой.  В случае кривой третьего порядка таких векторов четыре, при этом кривая начинается в точке ${{\vec{P}}_{0}}$, заканчивается - в ${{\vec{P}}_{3}}$, и имеет первую производную в начале и в конце траектории направленную в направлении точек ${{\vec{P}}_{1}}$ и ${{\vec{P}}_{2}}$ соответственно. 

Ломаная линия,  интерполируется, как показано на графике выше, кривой Безье с параметрами: ${{\vec{P}}_{0}}=(0,0,0)$,${{\vec{P}}_{1}}=(0,0,2.5)$,${{\vec{P}}_{2}}=(0,0.7,0.2)$, ${{\vec{P}}_{3}}=(0,0.7,2.7)$. Одна единица, как прежде, соответствует двум километрам. Сравнения сигналов от источников с гладким ({\bfseries bezier}) и ломаным ({\bfseries linear}) каналом приведены на двух графиках ниже. 
\ilus{zagagulja/E_bezier_zagagulja.eps}{Поле Е для источников ломаной и гладкой формы.}{ris_11}
\ilus{zagagulja/H_bezier_zagagulja.eps}{Поле H для источников ломаной и гладкой формы.}{ris_12}
Из приведённых графиков видно, что отличия небольшие и локализованы вблизи пиковых значений сигналов.
Заметим одну особенность, связанную с излучением такого канала со ступенькой. Если при интерпретации излучения молниевого разряда, имевшего канал изогнутой формы, использовать модель с прямолинейным источником, получится, что по каналу друг за другом проходят два импульса. Эта простая модельная ошибка показывает, что отказ от учета формы молниевого разряда в исследованиях излучения молний может привести к неправильным результатам.

\section{Оценка границ применимости дипольных представлений излучателя.}
\parsec Зададимся вопросом, насколько поля нашего источника отличаются от полей вертикального электрического точечного диполя, расположенного в начале координат. Аналитические выражения, связывающие спектр полей и источника, общеизвестны, в этой работе они были взяты из источника \cite{bib_13}. Зная эти выражения можно написать дифференциальные уравнения, связывающие сигналы - наблюдаемые поля и сигнал – источник.

Дифференциальное уравнение для поля E:
	\[\left[\partial^2_t + \left(\frac{c}{r} \frac{sin^2(\theta) - 2cos^2(\theta)}{sin^2(\theta)}
\right)\partial_t + \left(\frac{c^2}{r^2}\frac{sin^2(\theta) - 2cos^2(\theta)}{sin^2(\theta)}\right)\right]P_e(t) = -\frac{2\pi\epsilon_0c^2r}{sin^2(\theta)}E_z\left(t + r/c\right)\]
для поля H:
	\[\left[\partial^2_t + \left(\frac{c}{r}\right)\partial_t\right]P_e(t) =  \frac{2\pi\epsilon_0c^2Z_0r}{sin(\theta)}H_y\left(t+r/c\right)\]

Пусть наблюдатель располагается на прямой (x = 9 км, y=0). Будем поднимать точку наблюдения и для каждой точки вычислять эквивалентный дипольный момент по уравнениям, приведённым выше. Если бы источником был настоящий точечный диполь, все полученные сигналы совпадали бы. Но источник сложный, и сигналы, в общем, не совпадают. Параметры тока $\alpha$ и $\beta$ равны $10^4$ и $10^5$. Замедления нет ($\gamma = 0$). Угасания амплитуды волны нет ($\xi (l)= 1$). Источник – вертикальный канал длиной 3 км.
\ilus{img_Idl_dist9km.eps}{Эквивалентный ток точечного диполя.}{ris_13}
Для оценки недипольности использовалась L2  норма. Для каждой точки наблюдения вычисляются эквивалентные токи точечного диполя, полученные по полям E и H в этой точке и во всех, что ниже неё. Вычисляется средний сигнал. Вычисляется L2 разница между каждым сигналом и средним. Усреднённая разница даёт величину в процентах, приведённую в верхнем левом углу графика:
1.7\%  - 0 м, 2.8\%  - 1.9 км, 4.6\% – 3.9 км, 21\% – 5.8 км.
Стоит заметить, что при поднятиях таких, что $tan(\theta) = \sqrt{2},~~\theta \approx 54.7^\circ$, параметр $\alpha$ уравнения становится отрицательным. Это – неустойчивость уравнения, и она означает, что для слишком высоких точек наблюдения в общем не подобрать эффективный ток источника – точечного диполя.

\section{Заключение.}
\parsec Был разработан метод вычисления полей вблизи источника в задачах о излучении молниевого разряда. Полученная теория является достаточно общей и описывает многие модели источника, применяемые в публикациях другими авторами. Возможность комбинировать различные виды геометрии канала, формы волны тока, скорости запаздывания и ослабления делает полученные результаты применимыми к большому числу различных ситуаций. Используемый спектральный подход естественным образом обобщается до использования коэффициентов передач, получаемых решением одночастотных задач распространения.

	  Создана программа, реализующая теоретические изыскания. Язык программирования — C. Программа устроена таким образом, что формат входных данных удобен для подключения её в качестве модуля к проекту, нуждающемуся в реализуемых программой вычислениях. Программа способна эффективно использовать многопроцессорные системы, вычисляя значения спектра для разных частот параллельно на нескольких процессорных ядрах. Многопоточность реализована с помощью стандарта OpenMP.  Практической проблемой является отсутствие кроссплатформенности: программа разрабатывалась в среде ОС GNU/Linux, работоспособность в других UNIX-подобных ОС не проверялась, для сборки в среде Wndows определённо потребуется переработка.

	Для некоторых моделей источников были вычислены поля, произведены исследования влияния параметров модели на форму сигналов. В частности изучению подверглось влияние  уменьшения скорости движения токового импульса на форму, а также изменения сигналов при сглаживании кусочно-линейной траектории. Кроме того, была дана оценка отличию полей сложного источника от поля точечного диполя.

	Разработанные средства позволяют изучать особенности полей грозового разряда в связи с задачами сопоставления параметров источника с наблюдаемыми экспериментально полями ближней зоны.



\section{Литература.}
\begin{thebibliography}{20}
\bibitem{bib_1} Rakov V.A. , Uman M.A."Review and Evaluation of Lightning Return Stroke Models Including Some Aspects of Their Application", IEEE Trans. on EMC, vol. 40, No. 4, November 1998, part II, Special Issue on Lightning, pp. 403-426. 
\bibitem{bib_2} Bruce C. E. R. and Golde R. H., “The lightning discharge,” J. Inst. Elect. Eng.—Pt. 2, vol. 88, pp. 487–520, 1941.
\bibitem{bib_3}	M. N. Plooster, “Shock waves from line sources. Numerical solutions and experimental measurements,” Phys. Fluids, vol. 13, pp. 2665–2675, 1970.
\bibitem{bib_4}	S. I. Braginskii, “Theory of the development of a spark channel,” Sov. Phys. JETP, vol. 34, pp. 1068–1074, 1958 (Engl. transl.).
\bibitem{bib_5}	E. I. Dubovoy, M. S. Mikhailov, A. L. Ogonkov, and V. I. Pryazhinsky, “Measurement and numerical modeling of radio sounding reflection from a lightning channel,” J. Geophys. Res., vol. 100, pp. 1497–1502, 1995.
\bibitem{bib_6}	A. S. Podgorski and J. A. Landt, “Three dimensional time domain modeling of lightning,” IEEE Trans. Power Del., vol. PWRD-2, pp. 931–938, July 1987.
\bibitem{bib_7}	R. Moini, V. A. Rakov, M. A. Uman, and B. Kordi, “An antenna theory model for the lightning return stroke,” in Proc. 12th Int. Zurich Symp. Electromagnetic Compat., Zurich, Switzerland, Feb. 1997, pp. 149–152.
\bibitem{bib_8}	J. E. Borovsky, “An electrodynamic description of lightning return strokes and dart leaders: Guided wave propagation along conducting cylindrical channels,” J. Geophys. Res., vol. 100, pp. 2697–2726, 1995. 
\bibitem{bib_9}	G. H. Price and E. T. Pierce, “The modeling of channel current in the lightning return stroke,” Radio Sci., vol. 12, pp. 381–388, 1977.
\bibitem{bib_10}	 L. Baker, “Return-stroke transmission line model,” in Lightning Electromagnetics. New York: Hemisphere, 1990, pp. 63–74.
\bibitem{bib_11}	Y. T. Lin, M. A. Uman, J. A. Tiller, R. D. Brantley, W. H. Beasley, E. P. Krider, and C. D. Weidman, “Characterization of lightning return stroke electric and magnetic fields from simultaneous two-station measurements,” J. Geophys. Res., vol. 84, pp. 6307–6314, 1979.
\bibitem{bib_12}	C. A. Nucci, G. Diendorfer, M. A. Uman, F. Rachidi, M. Ianoz, and C. Mazzetti, “Lightning return stroke current models with specified channel-base current: A review and comparison,” J. Geophys. Res., vol. 95, pp. 20 395–20 408, 1990.
\bibitem{bib_13}	Г.И. Макаров, В.В. Новиков, С.Т. Рыбачек. Распространение электромагнитных волн над земной поверхностью. Москва ``Наука'' 1991. Страницы 28-29. 
\bibitem{bib_14}	Borisov V.V., Usupov I.E. Space structure of electromagnetic field components for different types of lightning channels. Proc. 24 Int. Conf. on Lightning Protection, Birmingham, 1998.
\end{thebibliography}

\section{Приложение 1.}
Применение дифференциальных операторов к потенциалу Герца
Связь между (11) и (14), (15) требует аккуратного вычисления. В выражениях, приведённых ниже, используется суммирование по повторяющимся индексам без знака суммы. Начнём с поля $E$. Выполним подстановку (11) в (12):
	\[E_{i\omega} = k^2 \frac{i}{4\pi\omega\varepsilon} \int\limits_{\sim\sim}dl ~~ \tau_i(l) I_w(l)\frac{e^{ikR}}{R} + 
\frac{i}{4\pi \omega\varepsilon}\int\limits_{\sim\sim}dl ~~ \tau_j(l)I_w(l)\partial_i\partial_j\frac{e^{ikR}}{R} =\]
	
	\[= \frac{ik}{4\pi\omega\varepsilon} \int\limits_{\sim\sim}dl ~~  I_w (\tau_i k^2 + \tau_j\partial_i\partial_j)\frac{e^{ikR}}{kR}  =\]
Воспользуемся выражением для второй производной $\frac{e^{ikr}}{kR}$, которое будет вычислено после:
	\[= \frac{ik}{4\pi\omega\varepsilon} \int\limits_{\sim\sim}dl ~~ I_w \left( \tau_i k^2 \frac{e^{ikR}}{kR} + \tau_j \left[ \frac{ie^{ikR}}{R^2} \left(\delta_{ij} - 3 n_i n_j \right)\left( 1 + \frac{i}{kR} \right) -  \frac{k e^{ikR}}{R} n_i n_j \right] \right) =\]
	\[= \frac{ik}{4\pi\varepsilon\omega} \int\limits_{\sim\sim}dl ~~ I_w(l) \left[ \tau_i k^2 \frac{e^{ikR}}{kR} + \frac{i e^{ikR}}{R^2} \left( \tau_i - 3 n_i(\vec{\tau},\vec{n}) \right)\left( 1 + \frac{i}{kR}\right) - \frac{k e^{ikR}}{R} n_i (\vec{\tau}, \vec{n}) \right]\]
Поле $E$ сосчитано.

Поле $H$, исходя из (11) и (13):
	\[H_{wi} = -i\omega ~ \varepsilon_{ikl} ~ \partial_k  \frac{i}{4\pi\omega}  \int\limits_{\sim\sim}dl ~~  dl I_w(l)\tau_l \frac{e^{ikR}}{R}\]
Здесь $\varepsilon_{ijk}$ - символ Леви-Чивиты. Он равен нулю, если есть повторяющиеся индексы. При перестановке двух стоящих вблизи символов его значение приобретает обратный знак. $\varepsilon_{123} = 1$.
	\[= \frac{1}{4\pi} \varepsilon_{ikl} ~ \int\limits_{\sim\sim}dl ~~ I_w(l)\tau_l \partial_k \frac{e^{ikR}}{R} =\]
	\[= \frac{1}{4\pi} \varepsilon_{ikl} ~ \int\limits_{\sim\sim}dl ~~ I_w(l)\tau_l \frac{e^{ikR}}{R}\left(ik - \frac{1}{R}\right) n_k =\]
Из символа Леви-Чивиты и двух векторов собирается векторное произведение:
	\[= \frac{1}{4\pi} \int\limits_{\sim\sim}dl ~~ I_w(l) \frac{e^{ikR}}{R}\left(ik - \frac{1}{R}\right) ~ [\vec{n},\vec{\tau}]_i\]
Поле $H$ сосчитано.
Вычисление производной $\partial_\eta\partial_\xi \frac{e^{ikR}}{kR}$.
Переменные индексы $\eta$, $\xi$ обозначают координаты $\{x,y,z\}$. Значения $R$, $n$ определяются выражениями (16).
	\[\partial_\eta\partial_\xi \frac{e^{ikR}}{kR} =\]
Пользуемся тем, что ${{\partial }_{\xi }}R={{n}_{\xi }}$:
	\[= \partial_\eta \frac{e^{ikR}}{kR}\left(ik - \frac{1}{R}\right) n_\xi =\]
Раскрываем $n_\xi = \frac{R_\xi}{R}$:
	\[= \partial_\eta \frac{e^{ikR}}{kR^2}\left(ik - \frac{1}{R}\right)R_\xi =\]
Обозначим часть выражения, не содержащую $R_\xi$, $f(R)$. По правилам дифференцирования, $\partial_\eta f(R)R_\xi = f'(R) n_\eta R_\xi + f(R)\delta_{\eta\xi}$. Применим это:
	\[= \frac{\partial}{\partial R} \left[ \frac{e^{ikR}}{kR^2}\left(ik - \frac{1}{R}\right) \right] n_\eta R_\xi +
 \frac{e^{ikR}}{kR^2}\left(ik - \frac{1}{R}\right)\delta_{\xi\eta}  =\]
	\[= \frac{k e^{ikR}}{R^2} \left( \frac{3}{k^2 R^2} - \frac{3i}{kR} - 1 \right) \frac{R_\eta R_\xi}{R} +
  \frac{e^{ikR}}{kR}\left(ik - \frac{1}{R}\right) \frac{1}{R} \delta_{\xi\eta} =\]
	\[= \frac{k e^{i k R}}{R} \left( \frac{3}{k^2R^2} - \frac{3i}{kR} - 1 \right) n_\eta n_\xi + 
	\frac{i e^{ikR}}{R^2} \left(1 + \frac{i}{kR}\right)\delta_{\xi\eta} =\]
	\[= \frac{3 k e^{i k R}}{R} \left( \frac{1}{k^2R^2} - \frac{i}{kR}\right) n_\eta n_\xi - 
\frac{k e^{i k R}}{R} n_\eta n_\xi + 
	\frac{i e^{ikR}}{R^2} \left(1 + \frac{i}{kR}\right)\delta_{\xi\eta} =\]
	\[= \frac{ - 3 i e^{i k R}}{R^2} \left( \frac{i}{kR} + 1\right) n_\eta n_\xi - 
\frac{k e^{i k R}}{R} n_\eta n_\xi + 
	\frac{i e^{ikR}}{R^2} \left(1 + \frac{i}{kR}\right)\delta_{\xi\eta} =\]
	\[= \frac{ie^{ikR}}{R^2} \left( \delta_{\xi\eta} - 3 n_\xi n_\eta\right)
	\left(1 + \frac{i}{kR}\right) - \frac{k e^{ikR}}{R} n_\xi n_\eta\]
В таком виде выражение готово к подстановке в формулы.
\end{document}
